Astronomy & Astrophysics manuscript no. 3673 


February 5, 2008 


(DOI: will be inserted by hand later) 





in 
o 
o 

(N 

Oh 

<D 



> 

o\ 
o 
in 
o 

S3 

Oh 
6 

03 



X 



Estimation and reduction of the uncertainties in chemical 
models: Application to hot core chemistry 

V. Wakelam 1 , F. Selsis 2 , E. Herbst 1 ' 3 , P. Caselli 4 

1 Department of Physics, The Ohio State University, Columbus, OH 43210, USA 

2 Centre de Recherche Astronomique de Lyon, Ecole Normale Superieure, 46 Allee d'ltalie, F-69364 Lyon cedex 
7, France 

3 Departments of Astronomy and Chemistry, The Ohio State University, Columbus, OH 43210, USA 

4 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi, 50125 Firenze, Italy 

Received 21 June 2005 / Accepted 25 August 2005 

Abstract. It is not common to consider the role of uncertainties in the rate coefficients used in interstellar gas- 
phase chemical models. In this paper, we report a new method to determine both the uncertainties in calculated 
molecular abundances and their sensitivities to underlying uncertainties in the kinetic data utilized. The method is 
used in hot core models to determine if previous analyses of the age and the applicable cosmic-ray ionization rate 
are valid. We conclude that for young hot cores (< 10 4 yr), the modeling uncertainties related to rate coefficients 
are reasonable so that comparisons with observations make sense. On the contrary, the modeling of older hot 
cores is characterized by strong uncertainties for some of the important species. In both cases, it is crucial to take 
into account these uncertainties to draw conclusions from the comparison of observations with chemical models. 
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1. Introduction 

Chemical models usually include thousands of reactions. 
The rate coefficients of these reactions are known to have 
uncertainties that affect the theoretical abundances com- 
puted by models. Estimating the resulting error bars is 
the only rigorous way to compare modeled and observed 
abundances. Although such an estimation of uncertainties 
has been undertaken in other fields where chemical net- 
works are used, e.g. i n atmospheric photochemistry (see 
iDobriievic et all2 003) , this important aspect has been rel- 
atively neglected in interstellar astrochemistry. To the best 
of our knowled ge, only thre e previous studies have ap- 
peared. In 1976 iKuntz et all published a sensitivity anal- 
ysis of molecular abundances to uncertainties in rate coef- 
ficients under steady-state conditions in dark clouds. The 
main purpose of this study was more the identification 
of chemical schemes forming the species than an estimate 
of "real" uncertainties in abundance. The two more re- 
cent works are by (i) iRoueff et al . (1996), who studied 
uncertainties in chemical modeling i n the region of bista- 
bility, and fii) IVasvunin et all {2004 h who also performed 
a study for the steady-state chemistry occurring in cold 
dark clouds. Independently we developed some procedures 
to include in a systematic manner the imprecision due 



to chemical data uncertainties in our interstellar chemi- 
cal models. In addition, a sensitivity analysis is utilized 
to focus on the major reactions that determine the uncer- 
tainties. In this paper, we study the implications on the 
modeling and dating of protostellar hot cores. 

During the collapse of a protostar, the temperature 
and density nearest the center increase to form a hot-core 
region where the dust mantles evaporate, releasing mole- 
cules into the gas phase. Considering the time of evapo- 
ration of the mantles as the initial time, hot cores can be 
dated by their chemical evolution. Indeed, the abundances 
of some chemical compounds evolve more rapidly than 
the structure of the protostellar envelope. Therefore, the 
comparison of the time-dependent abundances computed 
for these species with the observed ones can constrain the 
age of the source. Sulphur-bearing species are usually con- 
sidered to be good chemical clocks for hot cores: sulphur 
chemistry is initiated by evaporation from icy mantles and 
evolves sufficiently rapidly for the purpose jCharnlevll 997: 
lHatchell et al.lll998llWakelam et alJl20044 - The chemical 
network and model we use for this study has been pre- 
viously applied to sulphur chemistry in the hot corino 1 
of IRAS16293-2422 in order to constrain both the age of 
the source and the form of sulphur evaporated from the 
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grams <|Wakelam et all 120044 - The observed abundance 
ratios of sulphur-bearing species (SO2/SO, SO2/H2S and 
OCS/H2S) were reproduced, to within the observational 
error bars, only for an age of 1500—2500 yr, with the added 
assumption that the main source of evaporated sulphur is 
in the atomic form. We have now revisited these results by 
including uncertainties in the chemical rate coefficients. 

In this article, we first describe the general method to 
estimate the uncertainties in computed abundances and 
to identify the reactions of the chemical network that are 
mainly responsible for these uncertainties. Then, we show 
how the imprecision of the model results affects the scien- 
tific conclusions for the specific case of hot core chemistry. 
In particular, we explore the consequences in the use of 
molecular abundances to constrain the cosmic-ray-induced 
ionization rate and the age of a hot core. 
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Fig. 1. Distribution of the rate coefficient of the reaction 
O + SO — » SO2 + hf for 400 runs at a temperature 
of 100 K. The solid vertical line is the standard rate of 
UMIST. 



2. Chemical modeling 

2.1. Model Description 

We have use d the pseudo-time-depe ndent chemical model 
described in lWakelam et al.1 l|2004a|) . This code computes 
the chemical evolution of gas-phase species at a fixed gas 
temperature and density and for initial molecular abun- 
dances. The model includes 930 reactions involving 77 
species that contain the elements H, He, C, O and S. The 
most complex molecule in our reduced sample is proto- 
nated methanol (CHjOH.^), which contains seven atoms. 
The model includes 780 gas-phase reactions (ion-neutral, 
neutral-neutral and dissociative recombination processes), 
and allows species to stick to the surfaces of dust particles, 
as well as to evaporate and undergo charge exchange with 
the grains. The chemical network originally used for the 
mode l was taken from the osu.2003 network 2 l|Smith et alJ 
l2004h . For this study, however, we mainly used the rate 
coefficients given in the UMIST99 database 3 because it 
provides estimates for uncertainties in the reaction rate 
coefficients (see next section). This is the reason for small 
differences between th e calculated abundances given in 
I Wakelam et al.1 l|2004a|) and in the present study. We an- 
ticipate that these slight discrepancies remain within the 
estimated er ror at the considered te mperatures. 

Following lWakelam et al.l ( 2004af) , the initial gas phase 
composition is computed from the model for a molecular 
cloud with a temperature of 10 K and an H2 density of 
10 4 cm" 3 . The initial time (t=0) corresponds to the for- 
mation of the hot core, when the temperature increases 
to 100 K (the evaporation temperature of H 2 0), injecting 
into the gas phase the grain-mantle composition. We took 
the mantle compositio n observed in the env ironment of 
massive protostars Csee lWakelam et al"ll2004a|) . The abun- 
dance and form of sulphur initially evaporated from the 
grai ns are not quite known and remain a serio us problem 
(see iR.uffle et~al]ll999t IWakelam et all l2004ah . However, 

2 http:/ /www.physics. ohio-state.edu/~eric/research_files/ 



the species OCS and H2S seem to be present in the grain 
mantles. Indeed the first one has been observed in the solid 
state wit h a low fractional ab undance of 10 -7 compared 
with H 2 i Palumbo et alj fl997). Although the second one 
has never been detected on grains, the large abundances 
(> 10~ 8 with respect to H2) observed in hot cores and 
along outflows cannot be produced by gas phase routes 
(which accounts for less than 1% of the observed abun- 
dances). Thus, we assume that H2S is present on grains 
with an abu ndance lower than or equal t o 10 with re- 
spect to H2 (van Dishocck & Blake 1998). In addition to 
these molecules, Wakela m et al. l|2004a|) showed that the 
only possibility to reproduce the sulphur-bearing abun- 
dances observed in the low mass hot core of IRAS16293- 
2422 is to evaporate a large amount of atomic sulphur. 
Based on this prior analysis, we used a mantle composi- 
tion here in which H2S, OCS and S have abundances (with 
respect to H2) of 10~ 7 , 10 -7 , and 3 x 10~ 5 , respectively. 

2.2. Method 

To include uncertainties in kinetic data and to estimate 
their impact on the modeled abundances, we applied a 
Monte-Carl o procedure developed by Dobrijevic et al. 
<ll998_l2003h to study the gas-phase chemistry of the atmo- 
spheres of giant planets. This method was itself inspired by 
earlier studies dedicated to terrestrial stratospheric chem- 
istry dStolarski et alJ Il978t iThomoson fc Stewart! Il99lt 
IStewart fc Thompson! Il99ffll . First of all, we assume that 
the errors in the model results are due only to uncertain- 
ties in the gas-phase reaction rate coefficients, which im- 
plies that temperature, gas and grain densities, ionizing 
radiation, elemental abundances, and initial concentra- 
tions are all well characterized. This is obviously an ideal 
case but it allows us to evaluate the error specifically due 
to the kinetic data and to indentify the main reactions 
generating this error. 

Each reaction i included in the model possesses a rate 
coefficient K%. The UMIST kinetic database gives us the 
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direct cosmic-ray ionization, and cosmic-ray-induced 
photoreactions. The UMIST database also provides a 
factor Fi quantifying the uncertainty of the Ki. Under the 
physical conditions of the interstellar medium, this factor 
Fi itself cannot be estimated to a high level of precision, 
and rate coefficients in the UMIST database are classified 
in only 4 categories of precision: within 25% (F — 1.25), 
50% (F = 1.50), a factor of 2 (F = 2.0) and a factor of 10 
(F = 10). For the few reactions for which the coefficient 
F is not available, we assume an error within 25%. This 
choice is arbitrary but motivated by several arguments. 
To begin with, there is no obvious reason to assume that 
the reactions from osu.2003 but absent from UMIST are 
affected, on average, by a larger uncertainty than the 
others. Therefore, we chose the precision that dominates 
the set of UMIST reactions (40% of them have F = 1.25). 
Another reason has to do with the next goal of our study, 
which is to identify the main reactions responsible for the 
error in the modeling. By assuming a high uncertainty 
on the reactions without known F factor, our selection 
method would tend to point to these reactions. We prefer 
to avoid such bias, as we already know that these rate 
coefficients need to be considered with priority in order 
to have their uncertainty constrained. 

We assume (as illustrated in Fig. ^) that a rate coef- 
ficient K can be considered as a random variable, lognor- 
mally distributed over an uncertainty range centred on the 
recommended value. This basic ass umption is well ver ified 
in the absence of systematic error s ( Thompson fc Steward 
ll99ltlsTewart & Thompson! 19961) . For this reason, kinetic 
databases commonly provide the uncertainty as a A log K 
error (see for instance the IUPAC database: Atkinson et al. 
l2004h or as the factor F (which is equiva lent, as logF = 
A log K) . Here, our method differs from IVasvunin et alJ 
( 2004) since these authors used an equiprobable distribu- 
tion within the error interval given in UMIST. However, 
the authors noticed that their results are not significantly 
affected by the form of the distribution. 
Thus, given a factor Fi estimated at a la level of con- 
fidence, our approach implies that the value of log(i^i) 
follows a normal distribution with a standard deviation 
log(i 7 i). In other words, the probability to find Ki between 
K^i/Fi and K a<i x F; is 68.3%. 

With the standard set of rate coefficients Kqj, we can 
compute the standard evolution of the abundances Xoj(t) 
of the species j in the considered medium, which is what 
typical chemical models do. The specificity of our method 
consists of generating N new sets of rate coefficients (typ- 
ically N — 400) by taking into account the uncertainties 
affecting each Ki. This is done by generating N sets of 
random numbers e; with a normal distribution centered 
on and with a standard deviation of 1. For each set of 
e;, the new Ki are given by 



of the species j at time t. The mean value of logXj(f) 
(log Xj(t)) gives us the "recommended" value while the 
dispersion of log Xj(t) around logX, (i) determines the 
error due to kinetic data uncertainties. 



Once the error in the computed abundances is esti- 
mated, it is extremely useful to identify which reactions 
among the many possible are responsible for most of this 
error. It is indeed of high interest to be able to point out 
a few reactions for which more accurate measurements, or 
theoretical estimates, of the rate coefficient would signifi- 
cantly reduce the error in model results. 

Let us consider a species j at a time t. The error 
in its abundance due to uncertainties in kinetic data is 
AlogX,(t). Each of the Nr reactions included in the 
model contributes to this error. In order to estimate their 
individual contribution to the errors, we perform Nr new 
runs of the model. In each run i, we replace the standard 
rate Kqj of one reaction by its la upper value Ko,i x Fi, 
all the other rates being fixed to their standard values 
Kq. Each run i performed with such a perturbation of the 
rate Ki produces an abundance Xj(t) for every molecular 
species j. The ratio R l j(t), defined by the relation 



R)(t) = 



\X){t) 



(2) 



log Ki = log K 0A + e l log Fi 



(1) 



yields an index quantifying the influence of the reaction i 
on the total error for species j. The ratio Rj(t) can also 
be computed for abundance ratios instead of abundances 
when necessary. 

The rigorous way to estimate the individual error in- 
duced by a reaction i would be to compare the error 
A log Xj(t) obtained when all reaction rates are randomly 
chosen within their uncertainty range, with the one ob- 
tained when all the reactions but i are randomly cho- 
sen. The difference betwen these two values of A log Xj (i) 
would give the error in Xj(t) that is due to reaction i 
only. However, such an approach requires Nr x N runs 
(780 x 400 = 312,000 in our case), which is not afford- 
able in terms of computational time and data storage. 
Using the index Rj(t) is an alternative practical method 
although, in some cases, it may not be sufficient to iden- 
tify some of the reactions having a significant impact on 
the error. This is due to the non-linearity of chemical net- 
works and the complex structure of the routes to form 
some species. Nevertheless, once we have selected a group 
of reactions with the highest index Rj (t) we can artificially 
nullify or reduce their uncertainty and check the resulting 
decrease of the error AlogX, (t). It is important to note 
that a set of reactions identified to be responsible for most 
of the error in the abundance of one given compound j at 
a time t can have a minor effect on the error in the abun- 
dance of another species and/or at another time. 

Our approach c an be compared with the earlier work of 
iKuntz et alJ l|l976|) . who studied the sensitivity of molec- 
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Table 1. Conversion between A log A, relative abun- 
dances and error domains in X. 



3 4 5 

log[t(yr)] 



3 4 5 

log[t(yr)] 



Fig. 2. Abundances of SO and S computed with 400 runs 
(lower panel) as a function of time. The upper panels rep- 
resent histograms of the abundance distribution of these 
species for a given time (10 4 and 10 5 yr for SO and S 
respectively). The gas temperature is 100 K and the H2 
density 10 7 cm" 3 . 



31 reactions, allowing them to us e a Fourier analysis tech- 
nique (FAST, ICukier et alJll978h . which consists of peri- 
odically varying all the rate coefficients and solving the lin- 
ear system as a function of the coefficients. This method is 
efficient for small systems but the number of required solu- 
tions increases steeply with the number of reactions. With 
our network of 780 reactions, this technique would re- 
quire much more computation time than the Monte Carlo 
method we applied, without providing more information. 
Moreover, FAST is applicable only at steady state, which 
is not reached in astronomical objects such as dark clouds 
and hot cores, where time dependence has to be included. 
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probability of finding log A within log X±a and log A±2ct 
is 68.3% and 95.4%, respectively. 

Close to a strong inflection point in the evolution of 
the abundance (i.e., a steep edge in the time derivative 
of the abundance) , the distribution of log X is no longer 
Gaussian. Indeed, modifying the rate coefficients here pro- 
duces a temporal shift in the evolution of log X that can 
result in two separate distributions at the same time i, one 
on either side of the mean value log A. As an example, con- 
sider the inflection in the evolution of the atomic sulphur 
abundance in Fig.|5]and the resulting bimodal distribution 
at 10 5 yr. In such a case, the standard deviation seriously 
overestimates the error. Rigorously, when this occurs, we 
should fit the distribution by two Gaussians and define 
two disconnected error bars but, instead, we used the fol- 
lowing method to estimate AlogX at any time t. First, 
we calculate the normalized density of abundance- vs-time 
curves (or density of probablity) s ^ n x , where 5n is the 
number of curves per S log X interval and N is the to- 
tal number of runs. Then we identify the smallest interval 
[log A m i n , log A max ] (see Fig.^J) that contains 95.4% of the 
curves, and we define the error as follows: 



3. Results of Uncertainty Calculations 

In this section, we present general results concerning the 
uncertainties of computed abundances (with respect to 
H2) in hot cores for typical parameters: a temperature of 
100 K, a density of 10 7 cm" 3 , and ages lower than 10 6 yr. 
The two last parts of this section present the consequences 
of the rate coefficient uncertainties on two practical appli- 
cations: the determination of the H 2 cosmic-ray ionization 
rate and the determination of the hot core age using ob- 
served chemical abundances. 



3.1. Definition of the error 

The distribution of log X (the logarithm of the abundance 
with respect to H2 at an integration time t) is well fitted, 
at most integration times, by a Gaussian function centred 



AlogA = -(logA n 



logA min ) 



(3) 



When log A does follow a Gaussian distribution, this ex- 
pression reduces to A log A = 2a. In the rest of the 
paper, we refer to the error computed with this 
method, which has a level of confidence of 95.4%. 

Note again that abundance ratios can be treated the same 
way as abundances with respect to H2. 

If A log A = logs, the error domain [A m i n ,Amax] 
is [— , sA]. An alternative means of interpreting the er- 
ror is [A m ; n , Amaxl = A^a v + with A y~ = 1 — k and 



1. When A log A < 1. 



X 

AA_ 



on the mean value log A (see the abundance of SO in of X and [A m ; n , Jmax] = [ 1 023 1 1-023 A] . Tablefjgives the 



x — ^ ^-.^ 6 ^>. ^ -l, -~ AA + and 

the error A log A is proportional to the relative error in 

A : A log A w kHfix ' When A is lognormally distributed 

s = 10 2cr . As an example, if log A is known to within 

A log A±0.01, the error [AA_, AA+] in A is about ±2.3% 

x 
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Fig. 3. Evolution of the fractional abundances (/H2) of 
some species. The central solid line shows the evolution of 
the mean value logX while the 2 other lines delimit the 
domain containing 95.4% of the computed values of logX. 
The dashed lines give the 2ct standard deviation from the 
mean value with a Gaussian fit. Grey levels represent the 
density of probability (see text). 
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Fig. 4. Error AlogX of the abundance as a function of 
time for some species abundant in hot cores. 

3.2. Uncertainties in the abundances 

Figure |3| shows the evolution of logX computed for a va- 
riety of species at a temperature of 100 K and a H2 den- 
sity of 10 7 cm -3 . The evolution of logX ± AlogX (en- 
velopes containing 95.4% of the values as described above, 
solid lines on Fig. and \ogX ± 2a (standard devia- 
tion, dashed lines on Fig. |3J) are compared. Discrepancies 
between solid and dashed lines indicate a non-Gaussian 
distribution of logX. Figure 01 shows the evolution of 
the error A log X (estimated using the method explained 
above, equation^ as a function of time for some abundant 
species (X > lO^ 11 ) under the same physical conditions. 

The maximum A log X found in all our calculations 
is about 2 (for the species H2CO, CH3OH and S around 
10 5 yr). With AlogX = 2, the abundances X go from 
X/100 to lOQxX, which means that the total error spans 4 
orders of magnitude. However, such large error bars occur 
at late times when the hot core has probably ceased to 
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Fig. 5. The left panel shows the mean abundance of the 
SO2 and H2S molecules as a function of time (solid lines) 
with their error bars (grey contours) and the standard 
abundance (dashed grey lines). The right panel represents 
the ratio between the mean and the standard abundances 
(X/X ) for the two molecules as a function of time. 



the uncertainties on observed abundances produced solely 
by telescope and atmospheric calibrations. 

The error A log X in the computed fractional abun- 
dance depends on the species, as well as the integra- 
tion time (see Fig. 0J , the temperature, and the density. 
However, general patterns can be noticed. The error is log- 
ically related to the generational rank of the molecule: the 
abundances of the first generation of molecules at small 
integration times are determined by the initial conditions 
and are not affected seriously by kinetic data. For these 
species, the statistical error starts to be significant once 
the abundance, which typically decreases, is modified by 
the chemistry. This is the case for CH4, H2CO, CH3OH, 
H2S, OCS and S. Second or later generations of com- 
pounds are affected by the kinetic data uncertainties from 
the beginning, and the errors for their abundances gen- 
erally increase with increasing generational number. This 
effect can be explained by the increasing number of re- 
actions and species involved in their formation. The in- 
crease of the error with the generational rank can, how- 
ever, be verified only for species storing a small fraction 
of the global reservoir of the chemical elements of which 
they consist. Indeed, because the elemental abundances 
in the gas phase are fixed, the sum of the abundances 
of all species bearing a given element represents a con- 
stant reservoir that is not affected by uncertainties due 
to kinetic data. For instance, A [Sulphur] =0 , where the 
elemental abundance of sulphur is given by 



[Sulphur] = [S] + [SO] + [S0 2 ] + [CS] + [HS] + [H 2 S] - 



•(4) 



For t < 10 yr, atomic S is the dominant S-bearing species 
by two orders of magnitude, and the error of its abundance 
is thus negligible. At t > 10 5 yr, most of the sulphur in 
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Fig. 6. Abundance ratios (solid central lines) with er- 
ror bars computed using two methods: A\og(A/ B) (solid 
lines) and AlogA + AlogB (dashed lines). 



comes insignificant. In their study of error propagation 
in the chemical mode ling of interstellar molecular clouds, 
IVasvunin et alJ (|2004^> found a correlation between the er- 
ror and the number of atoms per species. We did not find 
such correlation in the hot corino chemistry. This might 
be due to differences in the physical conditions but also to 
the fact that we used a reduced sample of species. Ninety 
percent of the species included in our model have four 
atoms or fewer, and only four molecules possess six or 
seven atoms. Therefore, we cannot determine if such a 
correlation exists for the complex species with seven or 
more atoms, although it seems likely since these species 
are expected to have low abundances and be of late gen- 
erational rank. 

For the same physical conditions, the average molecu- 
lar abundance, X, is normally close to the standard abun- 
dance X (computed with the standard rates). However, 
for some species, a factor of up to 2 difference is found 
between the two values, especially when the distribution 
of X diverges from a log-normal distribution. This is for 
instance the case for H 2 S (as seen on Fig. EJ) around 
t=6 x 10 4 yr. On the same figure, we can see that the 
fast decrease of the H 2 S abundance produces sharp spikes 
in the ratio X /Xq, which occur when the curves of X and 
Xo cross and re-cross each other. 

3.3. Uncertainties in abundance ratios 

The uncertainty in abundance ratios can be computed us- 
ing two different approaches. The first one involves the cal- 
culation of the abundance ratios for each run and the de- 
termination of the error A\og(A/ B). The second method 
involves the error for each abundance and their subsequent 
addition: A [log A] + A [log B] . Fig.[|)]shows the evolution of 
the logarithm of some interesting abundance ratios with 



similar. This is the case for SO2 and H2S because, in our 
model, S0 2 is formed by the oxygenation of the initial 
atomic sulphur more than from the destruction of H2S. 
Another example is the ratio H2CS/OCS before 10 4 yr 
since OCS only starts to be destroyed after this time to 
form H2CS. When the species A and B are chemically re- 
lated, A [log A] + A[log_B] overestimates the error on the 
ratio. The abundance ratios CS/SO and SO2/SO are two 
good examples since the formation mechanism of CS and 
S0 2 are directly related to SO. For H 2 CS/OCS the dis- 
crepancy becomes huge for integration times between 10 5 
and 10 6 yr: A [log A] + A[logB] covers 7 orders of magni- 
tude while A[log(A/_B)] remains between 2 orders of mag- 
nitude. As a consequence in this paper, we use the error 
computed for the abundance ratios A[log(A/£?)]. 



4. Some consequences for hot core models 

4.1. The cosmic-ray ionization rate in IRAS 16293-2422 

The ionization rate £ is an important parameter of chem- 
ical models but is weakly constrained. It depends on the 
cosmic-ray flux, which cannot be directly determined. 
Some constraints have been derived from the comparison 
of observed abundance s with numerical predictions done 
with s everal values of C t f Wootten et all 19 791: C^gglh^g^l 



1998; Ivan der Tak fc van Dishoeckl 2000: ICaselli et al 
l2002HDotv et alJl2004l) . None of the models used included 
any uncertainties in the rate coefficients. Here we address 
the question whether or not it is possible to distinguish 
among several ionization rates in hot cores. 

For that purpose, we ran models for hot cores (100 K 
and 10 7 cm -3 ) with three different values of the ionization 
rate: ( = 1.3 x 10" 17 , 5.3 x 10" 17 and 1.3 x 10" 16 s~ x . 
Then, we searched for those species with abundance or 
abundance ratio that would allow us to distinguish among 
the values of £ despite the dispersion due to kinetic uncer- 
tainties. Our results show that the best molecules to con- 
strain the ionization rate in hot cores appear to be OH, 
CS and H3O 4 " (see Fig.[7|). Indeed, the abundance of OH 
remains almost constant and the abundances of CS and 
H30 + do not vary much with time before 10 4 yr. While 
OH can be detected at a large number of frequencies un- 
der 100 GHz thanks to its lambda-doubling transitions, 
H30 + can be detected via rotational- inversion transitions 
at significantly higher frequencies. If the age of the source 
is known, the OCS, H 2 S, and SO molecules can also be 
used. On the contrary, the species H2CS, CH3OH, H2CO, 
O2 and SO2 show distinct but close abundance distribu- 
tions only for the most extreme rates 1.3 x 10~ 17 and 
1.3 x 10 -16 s" 1 and intermediate values of the rate can not 
be distinguished. Other species, like H 2 0, CH4, C2H2, CO 
and O are less sensitive to ( at young ages (t < 10 5 yrs) 
and exhibit strong overlaps after. 

From both the theoretical and observational points of 
view, the use of abundance ratios instead of individual 



Wakelam et al.: Chemical modeling uncertainties 



7 




2 3 4 5 6 2 3 4 5 62 3 4 5 6 

log (t) log (t) log (t) 

Fig. 7. The upper panel presents calculated abundances (with respect to H2) of OH, CS and H30 + as a function of 
time with an error bar (95.4% confidence level) for three different ionization rates. The lower panel presents three 
calculated abundance ratios (H2CS/H2S, OCS/SO and H2CS/OCS) with similarly estimated uncertainties for the 
three cases. The solid lines in the lower panel represent the ratios observed in IRAS16293-2422 taking into account 
the uncertainties in the observations. 



on the H2 column density and the size of the source, as- 
suming that both molecules of the ratio come from the 
same region. Also, theoretical abundance ratios are less 
sensitive to the initial conditions, which are usually not 
well constrained. Our results show that some of the com- 
monly used abundance ratios involving sulphur-bearing 
species can then be used without confusion to distinguish 
among the three cosmic-ray ionization rates. As an ex- 
ample, we show the evolution of three of the abundance 
ratios: H 2 CS/H 2 S, OCS/SO and H 2 CS/OCS for the three 
values of C in Fig- 13 

The abundance of OH has never been determined in 
any hot core, to the best of our knowledge. Moreover, 
there exist only observations of low-energy transitions of 
CS (E up < 45 cm^ 1 ), which do not allow determination of 
its abundance in the inner regions of protostellar objects 
ijSchoier et all [2002). CS transitions with J > 7 cannot 
be observed with current ground based telescopes because 
of the receivers frequency ranges or they are absorbed by 
the atmosphere. There was one attempt to dete ct H3O 4 " 
towards IRAS16293-2422 bv lPhillips et all |l992) with no 
success. The authors deduced an upper limit of 2 x 10 -10 
for the abundance of this molecule which would strongly 
constrain the ionization rate in this source. However, this 
limit is for a 18" beam size. If we consider a 2" source 
size for the hot core of IRAS16 293-2422 and a H 2 co lumn 
density of 7.5 x 10 22 cm -2 (see lCeccarelli et al.ll2000h . the 
new limit on the HsO + abundance (compared with H 2 ), 
assuming LTE, is 4 x 10" 8 . The limit is thus too high to 
conclude anything when comparing with Fig. [3 Note that 
in hot core s, the typical "ionization tracers" (i.e. HC O + 
and CCH, lYan fc Dalgarnolll997t iGerin et al.lll997|) are 
not abundant. Our attention thus shifts to the ratios of 
the sulphur-bearing species discussed above. 



ties for the low mass protostellar source IR AS16293-2422 
l|Schoier et al.ll2002l:IWakelam et al.ll2004bh . The observed 
OCS/SO ratio clearly constrains the cosmic-ray ioniza- 
tion rate to be 1.3 x 10~ 17 s _1 . The two other observed 
ratios H 2 CS/H 2 S and H2CS/OCS are in agreement with 
this rate for ag e s arou nd 10 3 yr. As already noticed by 
Wa kelam et alJ ll2004afl . this result is contradictory to 
who could reproduce the abundances 



Dotv et all 



towards this source only with a hi gh cosmic-ray io nization 
rate of 1.3 x 10~ 16 s" 1 fsee lWakelam et al.ll2004al for dis- 
cussion) . Note that both the rate and the age were already 
determined bv l Wakelam et alJ l)2004al) and the goal here is 
to show that even with the introduction of the uncertain- 
ties in the reaction rates, we can still distinguish ionization 
rates between 1.3 x 10~ 17 and 1.3 x 10~ 16 s _1 . In other 
words, molecular abundances can be used to constrain the 
ionization rate, but one needs to verify that uncertainties 
in the reaction rate coefficients do not confuse the results. 



4.2. Consequence on the age of IRAS16293-2422 hot 
corino 



IWakelam et alJ l)2004af) constrained the age of the 
IRAS 16293-2422 hot corino by comparing the S0 2 /SO, 
S0 2 /H2S and OCS/H2S abundance ratios observed to- 
wards the source with the predictions of their chemical 
model. In their comparison, they only took into account 
the errors in the observed ratios , which are due to the 
calibration error of the telecopes. I Wakelam et alJ 1 2004al) 
obtained an age of (2±0.5) x 1 3 yr, where the uncert ainty 
in the age, 25% (see Fig. 7 of lWakelam et alJl2004a|) . We 
did the same comparison including the uncertainty in the 
rate coefficients with the method described in § 12.21 In 
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Fig. 8. Comparison between the observed abundance ratios SO2/SO, SO2/H2S and OCS/H2S (black empty boxes) 
towards IRAS16293-2422 and the theoretical predictions (gray contours). A: including the actual uncertainties in rate 
coefficients; B: decreasing the uncertainties of the 5 reactions in Tabled to 0; C: decreasing these uncertainties to 
achievable values (see Tabic 



Table 2. List of important reactions. 



Reaction 





P 


7 


UMIST uncertainties 


Achievable uncertainties 


O + SO -> S0 2 + hu 


3.20e-16 


-1.50 


0.0 


Error within factor 2 


20% 


H 2 + CR* -» H+ + e" 


1.30e-17 


0.00 


0.0 


Error within factor 2*** 


20% 


H 2 + CRPHOT** -> OH + H 


1.30e-17 


0.00 


971.0 


Error within factor 2 


20% 


CO + S -> OCS + hv 


1.60e-17 


-1.50 


0.0 


Error < 25 % 


25% 


H+ + S -> HS+ + H 2 


2.60e-09 


0.00 


0.0 


Error within factor 2 


10% 



*Cosmic-ray ionization; ** Cosmic-ray-induced photodissociation 

'"Uncertainty based on laboratory experiments rather than uncertainty in the cosmic-ray ionization rate 



modeling. In this case, we obtain an age between 800 and 
2, 500 yr, (1.65 ± 0.85) x 10 3 yr, i mplying an error of 50% 
in the age 4 . As already noticed in IWakelam etall l(2004al) . 
this age is relati vely short compared with previous esti- 
mates (~ 10 4 vrlCeccarelli et al-lEoOot iMaret et"al]l2002t 
Wakelam et al . 2004b). However, this determined age rep- 
resents only the time from the evaporation of the man- 
tles; the dynamical time needed to reach the required lu- 
minosity to form the observed hot core is an additional 
- 3 x 10 4 yr. 

Using the sensitivity analysis described in § 12.21 we 
then determined the reactions that produce most of the 
uncertainties in the three ratios (SO2/SO, SO2/H2S and 
OCS/H2S) at an age of 10 3 - 10 4 yr. We found five im- 
portant reactions, which are summarized in Table with 
their rate coefficients (a. (3, 7 from the UMIST database) 
and their uncertainties. The reactions are given in decreas- 
ing order of importance. Note that reactions 2 and 3 are 
directly related to the cosmic-ray ionization rate studied 
in § 14.11 and their uncertainties are based at least par- 
tially on astronomical rather than laboratory considera- 
tions. To check the importance of these five reactions, we 
decreased their uncertainties to and to achievable val- 
ues (see Table |5J| as shown on Fig. [S] (panels B and C) . 
In estimating achievable values, we were helped by Nigel 



Table 3. Range of computed abundance ratios and age of 
IRAS16293-2422 hot corino at t=10 4 yr. 



4 The age and its relevant error are denned as follows: let us 
define At = [t min , t max ] as the interval of time where the three 
observed ratios are reproduced within the uncertainties by the 





[minimum 


value, maximum value] 




current 


null 


achievable 


S02/SO 


[1.12, 2.34] 


[1.60, 1.64] 


[1.45, 1.82] 


S02/H2S 


[25.1, 47.9] 


[28.8, 41.7] 


[28.8, 41.7] 


OCS/H2S 


[11.2, 17.0] 


[12.0, 15.8] 


[11.7, 16.2] 


Age (yr) 


[800, 2500] 


[890, 1770] 


[890, 1990] 



Adams (private communication) . The uncertainties for the 
radiative association reactions between neutral species are 
rather speculative and contain the assumption that exper- 
iments to determine the rate coefficients can have the same 
low uncertainty as other experiments in neutral-neutral 
reactions. The uncertainty in panels B and C of Fig. [S] 
is still high at later times (~ 10 6 yr) because other re- 
actions than the ones listed in Table [21 are important. 
Table |21 gives the range of the computed abundance ratios 
SO2/SO, SO2/H2S, OCS/H2S at 10 4 yr using the current 
(col. 1), null (col. 2) and achievable (col. 3) uncertainties 
for the five reactions listed in Table [2 Note that in the 
case of these abundance ratios, we found that the relative 
error is well determined by the 2er dispersion (see § I3.1|l . 
In the last line of the table, we report the range of age 
found for the IRAS16293-2422 hot core. 



5. Conclusions 
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and some methods to improve the accuracy of theoretical 
predictions. We first present a general method to include 
the treatment of rate imprecisions in models and a pro- 
cedure to identify the main uncertain reactions responsi- 
ble for the error on the abundances. We found that for 
hot core models the relative errors A ^ on the predicted 
abundances of the more abundant species are lower than 
50% (with a 95.4% level of confidence) for times earlier 
than 10 4 yr, which is comparable to the uncertainties in 
observed abundances. This error depends on time, H 2 den- 
sity, and temperature, and is related to the generational 
rank of the molecule in the hot core chemistry. The first 
generation of molecules, which are evaporated from the 
grains, has lower uncertainties than the second and later 
generations. 

We also studied the consequences of rate uncertainties 
on the use of models to constrain the ionization rate £ 
and the age of hot cores. Among all the molecules, we 
found that OH, CS and fi30 + were the best to distinguish 
among £ = 1.3 X 1CT 17 , 5.3 x KT 17 and 1.3 x 1(T 16 s" 1 . 
Indeed, these molecules have reasonable abundances (> 
10 -11 ) which do not vary much with time before 10 4 yr 
and depend on £. Sulphur bearing abundance ratios such 
as H 2 CS/H 2 S, OCS/SO and H 2 CS/OCS can also be used 
since these molecules are easily observed in hot cores and 
show strong dependence on (. 

We compared our modeling including the uncertain- 
ties with observed S-bearing abundance r atios in order to 
constra in the age of IRAS16293-2422, as IWakelam et alJ 
ll2004ah did without the rate uncertainties. We also devel- 
oped a simple method to sort the reactions by the uncer- 
tainty they produce in the model results. This allowed us 
to identify a group of 5 reactions mainly responsible for 
the error in the S0 2 /SO, OCS/H 2 S and S0 2 /H 2 S ratios, 
at integration times between 10 3 and 10 4 yr. The decrease 
of the actual error for these 5 reactions to achievable val- 
ues would decrease the uncertainty of the age of this hot 
core from ±50% to ±38%. 

In conclusion, for young hot cores (< 10 4 yr), the mod- 
eling uncertainties related to rate coefficients are reason- 
able and comparisons with observations make sense. On 
the contrary, the modeling of older hot cores is charac- 
terized by strong uncertainties for some of the important 
species. In both cases, it is crucial to take into account 
these imprecisions to draw conclusions from the compari- 
son of observations with chemical models. In addition, be- 
ing able to identify, among the thousands of reactions in- 
volved in interstellar chemical networks, the few reactions 
for which a high accuracy is required can be useful espe- 
cially for old regions. Studies such ours rely on the esti- 
mation of the uncertainties provided by kinetic databases 
and thus, an effort should be done to better quantify and 
the rate coefficient uncertainties. 

Acknowledgements. V. W. and F.S. wish to thank Michel 
Dobrijevic for fruitful discussions on error propagation in 
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